Results (Julia)

Code
using RCall
include("main.jl")
WARNING: replacing module DiscontinuityIO.
L"Radial Distance ($AU$)"
Code
j_events_taus = load_taus(60:-10:20)

# find the numbers of events for each tau
combine(groupby(j_events_taus, :tau), nrow)
Number of events: 18088
Number of events: 20129
Number of events: 22606
Number of events: 26333
Number of events: 31793
5×2 DataFrame
Row tau nrow
Int64 Int64
1 20 31793
2 30 21007
3 40 14952
4 50 11904
5 60 9023
Code
w_events = load("data/events.Wind.fit.ts_0.09s_tau_60s.arrow");
┌ Warning: automatically converting Arrow.Timestamp with precision = MICROSECOND to `DateTime` which only supports millisecond precision; conversion may be lossy; to avoid converting, pass `Arrow.Table(source; convert=false)
└ @ Arrow /Users/zijin/.julia/packages/Arrow/5pHqZ/src/eltypes.jl:273

Check the discontinuities properties with time and radial distance

Code
datalimits_f = x -> quantile(x, [0.02, 0.98])
#32 (generic function with 1 method)
Code
df = w_events
plot_dist(
    data(df) * mapping(color = :year);
    datalimits=datalimits_f
)
easy_save("wind_distribution_time")
┌ Info: Saved /Users/zijin/projects/ids_spatial_evolution_juno/figures/wind_distribution_time.png
└ @ beforerr /Users/zijin/.julia/dev/beforerr.jl/src/utils/makie.jl:34
Code
plot_dist(
    data(df) * mapping(color = :year);
    maps = (density_log_map, B_log_map, b_log_map,),
    datalimits=datalimits_f,
)
easy_save("sw_paramters_time")
┌ Info: Saved /Users/zijin/projects/ids_spatial_evolution_juno/figures/sw_paramters_time.png
└ @ beforerr /Users/zijin/.julia/dev/beforerr.jl/src/utils/makie.jl:34
Code
df = w_events

plot_dist(
    data(df) * mapping(color = :year);
    datalimits=datalimits_f,
)
easy_save("wind_distribution_time")
Code
plot_dist(
    data(j_events_taus) * mapping(color=:r);
    datalimits=datalimits_f,
)
easy_save("juno_distribution_r")
┌ Info: Saved /Users/zijin/projects/ids_spatial_evolution_juno/figures/juno_distribution_r.png
└ @ beforerr /Users/zijin/.julia/dev/beforerr.jl/src/utils/makie.jl:34

The effect of tau

Code
plot_dist(
    data(j_events_taus_unique) * mapping(color=:label, row = :r);
    maps=(l_log_map, l_norm_log_map, j_log_map, j_norm_log_map, b_log_map,),
    datalimits=datalimits_f,
    figure_kwargs=(size=(1200, 1000),)
)

easy_save("juno_distribution_r_taus")
┌ Info: Saved /Users/zijin/projects/ids_spatial_evolution_juno/notebooks/figures/juno_distribution_r_taus.png
└ @ beforerr /Users/zijin/.julia/dev/beforerr.jl/src/utils/makie.jl:34

Archived (Following codes may not work)

Code
j_events_low_fit = load()
j_events_tau_20_fit = load(tau = 20)
j_events_high_fit = load(ts = 0.12)

# add a label column to the dataframes
j_events_low_fit.label .= "1 Hz (fit)"
j_events_high_fit.label .= "8 Hz (fit)"
j_events_tau_20_fit.label .= "1 Hz, 20s (fit)"
j_events = reduce(vcat, [j_events_low_fit, j_events_high_fit, j_events_tau_20_fit]) |> process;
Code
# NOTE: log axis for density is not working now
function plot_l_j_r(df)
    data_layer = data(df) * mapping(row=:r, color=:label)
    plot_dist(data_layer; datalimits=datalimits_f, figure_kwargs=(size=(1200, 1200),))
end
plot_l_j_r (generic function with 1 method)
Code
plot_l_j_r(j_events)
Code
function plot_l_j_r_hist(
    df;
    maps=(l_map, l_norm_map, j_map, j_norm_map),
    f_maps=(row=:r, color=:label),
    normalization=:pdf,
    fig_options = (size = (1200, 1000),)
)
    fig = Figure(;fig_options)

    plt = data(df) * mapping(; f_maps...) * (visual(Lines) + visual(Scatter))
    plt *= histogram(datalimits=datalimits_f, normalization=normalization)

    plts = [plt * mapping(m) for m in maps]
    axis = (; yscale=log10)

    grids = [draw!(fig[1:5, i], p, axis=axis) for (i, p) in enumerate(plts)]
    legend!(fig[0, 1:end], grids[1], titleposition=:left, orientation=:horizontal)

    fig
end
plot_l_j_r_hist (generic function with 1 method)
Code
plot_l_j_r_hist(j_events)
Figure 1
Code
maps = (l_log_map, l_norm_log_map, j_log_map, j_norm_log_map)
plot_l_j_r_hist(j_events; maps = maps, normalization = :none)
Code
function plot_l_j_r_hist(
    df;
    maps=(l_map, l_norm_map, j_map, j_norm_map),
    f_maps = (color=:r,),
    size = (1200, 800),
)
    fig = Figure(size=size)
    axis = (; yscale=log10)

    plt = data(df) * mapping(; f_maps...) * (visual(Lines) + visual(Scatter))
    plt_h = plt * histogram(datalimits=datalimits_f, normalization=:pdf)
    plts = [plt_h * mapping(m) for m in maps]
    grids = [draw!(fig[1, i], p, axis=axis) for (i, p) in enumerate(plts)]


    plt_h = plt * histogram(datalimits=datalimits_f)
    plts = [plt_h * mapping(m) for m in maps]
    grids = [draw!(fig[2, i], p, axis=axis) for (i, p) in enumerate(plts)]
    legend!(fig[0, 1:end], grids[1], titleposition=:left, orientation=:horizontal)

    fig
end

plot_l_j_r_hist(j_events_low_fit; maps=maps, f_maps=f_maps)
Code
@rput j_events
@rput j_events_taus
@rput j_events_taus_u
@rput j_events_low_fit
@rput j_events_tau_20_fit


R"""
library(ggplot2)
library(ggpubr)
library(patchwork)
source('utils.R')
"""
RObject{VecSxp}
$value
function (df1, df2, p1title = p1title, y_lim_1 = c(500, 40000), 
    y_lim_2 = c(1, 200)) 
{
    x_col <- "radial_distance"
    x_bins <- 5
    y_bins <- 12
    add_mode <- TRUE
    y_col <- "L_k"
    y_lim <- y_lim_1
    ylab <- lab_l_log
    p <- plot_binned_data(df1, x_col = x_col, y_col = y_col, 
        x_bins = x_bins, y_bins = y_bins, y_lim = y_lim, y_log = TRUE, 
        add_mode = add_mode)
    p1 <- p + labs(x = NULL, y = ylab) + ggtitle(p1title)
    y_col <- "L_k_norm"
    y_lim <- y_lim_2
    ylab <- lab_l_norm_log
    p <- plot_binned_data(df1, x_col = x_col, y_col = y_col, 
        x_bins = x_bins, y_bins = y_bins, y_lim = y_lim, y_log = TRUE, 
        add_mode = add_mode)
    p2 <- p + labs(x = x_lab_r, y = ylab)
    x_col <- "time"
    y_col <- "L_k"
    y_lim <- y_lim_1
    p <- plot_binned_data(df2, x_col = x_col, y_col = y_col, 
        x_bins = x_bins, y_bins = y_bins, y_lim = y_lim, y_log = TRUE, 
        add_mode = add_mode)
    p3 <- p + labs(x = NULL, y = NULL) + ggtitle(p2title)
    y_col <- "L_k_norm"
    y_lim <- y_lim_2
    p <- plot_binned_data(df2, x_col = x_col, y_col = y_col, 
        x_bins = x_bins, y_bins = y_bins, y_lim = y_lim, y_log = TRUE, 
        add_mode = add_mode)
    p4 <- p + labs(x = x_lab_t, y = NULL)
    (p1 + p2 + p3 + p4) + layout & scale_fill_viridis_c(limits = c(0.01, 
        0.25), trans = "log10", name = "pdf")
}

$visible
[1] FALSE
Code
R"""
density_plot <- function(df,
  x,
  density.y = "density",
  facet.by = "r",
  color = "label",
  strip.position = "top")
{
  p <- ggdensity(
    df, 
    x = x, y = density.y,
    color = color,
    facet.by = facet.by,
    add = "median",
    alpha = 0
  ) %>%
  facet(facet.by, ncol = 1, strip.position = strip.position)
}

plot <- function(
  df,
  d1_ulim = 8000,
  d2_ulim = 40,
  d3_ulim = 8,
  d4_ulim = 1.2,
  xscale = "none",
  yscale = "log10",
  p1_xlim = NULL,
  p1_ylim = NULL,
  p2_xlim = NULL,
  p2_ylim = NULL, # c(10^-3, 1)
  p3_xlim = NULL, p3_ylim = NULL,
  p4_xlim = NULL, p4_ylim = NULL,
  plot_func = density_plot,
  ...
) {

  p1 <- plot_func(
    filter(df, L_k < d1_ulim), "L_k", ...
  )
  
  p2 <- plot_func(
    filter(df, L_k_norm < d2_ulim), "L_k_norm", ...
  )
  
  p3 <- plot_func(
    filter(df, j0_k < d3_ulim), "j0_k", ...
  )
  
  p4 <- plot_func(
    filter(df, j0_k_norm < d4_ulim), "j0_k_norm", strip.position = "right", ...
  )

  p1 <- ggpar(p1,
      xscale=xscale, xlim = p1_xlim,
      yscale=yscale, ylim = p1_ylim
    ) + labs(x = lab_l)
  
  p2 <- ggpar(p2, 
      xscale=xscale, xlim = p2_xlim,
      yscale=yscale, ylim = p2_ylim, ylab=FALSE
    ) + labs(x = lab_l_norm)
  
  p3 <- ggpar(p3,
      xscale=xscale, xlim = p3_xlim,
      yscale=yscale, ylim = p3_ylim, ylab=FALSE
    ) + labs(x = lab_j)
  
  p4 <- ggpar(p4,
      xscale=xscale, xlim = p4_xlim,
      yscale=yscale, ylim = p4_ylim, ylab=FALSE
    ) + labs(x = lab_j_norm)
  
  tag_pool <- unique(df$r)
  
  p1 <- tag_facet(p1, open = "(a.", tag_pool = tag_pool)
  p2 <- tag_facet(p2, open = "(b.", tag_pool = tag_pool)
  p3 <- tag_facet(p3, open = "(c.", tag_pool = tag_pool)
  p4 <- tag_facet(p4, open = "(d.", tag_pool = tag_pool) + theme(
    strip.background = element_blank(),
    strip.text.x = element_blank()
  )

  p <- p1 + p2 + p3 + p4
  p <- p + plot_layout(nrow = 1, guides = "collect") & theme(legend.position='top')  
}
"""

Properties over radial distance with different tau values

Code
R"""
plot(
    j_events_taus,
    d1_ulim = Inf, p1_xlim = c(0,1e4), p1_ylim = c(1e-5, 1e-3),
    d2_ulim = 25, p2_xlim = c(0,25), p2_ylim = c(1e-2, 1e0),
    d3_ulim = 20, p3_xlim = c(0,8), p3_ylim = c(1e-2, 1e0),
    d4_ulim = 2, p4_xlim = c(0, 1), p4_ylim =  c(1e-1, 1e1),
)

save_plot("l_j_fit_taus", width = 15, height = 5)


plot(
    j_events_taus,
    xscale = "log10",
    d1_ulim = Inf, p1_xlim = c(1e2, 1e5), p1_ylim = c(1e-2, 1e0),
    d2_ulim = Inf, p2_xlim = c(1e-1, 1e2), p2_ylim = c(1e-2, 1e0),
    d3_ulim = Inf, p3_xlim = c(1e-2, 5e1), p3_ylim = c(1e-2, 1e0),
    d4_ulim = Inf, p4_xlim = c(1e-3, 5e0), p4_ylim = c(1e-2, 1e0),
)

save_plot("l_j_fit_taus_xscale_log", width = 15, height = 5)
"""
Code
R"""
plot(
    j_events_taus_u,
    xscale = "log10",
    d1_ulim = Inf, p1_xlim = c(1e2, 1e5), p1_ylim = c(1e-2, 1e0),
    d2_ulim = Inf, p2_xlim = c(1e-1, 1e2), p2_ylim = c(1e-2, 1e0),
    d3_ulim = Inf, p3_xlim = c(1e-2, 5e1), p3_ylim = c(1e-2, 1e0),
    d4_ulim = Inf, p4_xlim = c(1e-3, 5e0), p4_ylim = c(1e-2, 1e0),
    facet.by = NULL, color = "r"
)

save_plot("l_j_fit_xscale_log", width = 15, height = 5)

plot(
    j_events_taus_u,
    d1_ulim = Inf, p1_xlim = c(0, 1e4), p1_ylim = c(1e-5, 1e-3),
    d2_ulim = Inf, p2_xlim = c(0, 25), p2_ylim = c(1e-2, 1e0),
    d3_ulim = 20, p3_xlim = c(0,8), p3_ylim = c(1e-2, 1e0),
    d4_ulim = 2, p4_xlim = c(0, 1), p4_ylim =  c(1e-1, 1e1),
    facet.by = NULL, color = "r",
)

save_plot("l_j_fit", width = 15, height = 5)
"""

l_j_fit

l_j_fit_xscale_log

l_j_fit_xscale_log_count

l_j_fit_taus

l_j_fit_taus_xscale_log

Check the effect of time resolution and time window

Code
R"""
# plot(j_events_der)
# save_plot("l_j_r_der", width = 15, height = 10)
# plot(
#     j_events,
#     p1_ylim = c(1e-5, 1e-3),
# )
# save_plot("l_j_r_fit", width = 15, height = 10)

plot(
    j_events,
    xscale = "log10",
    d1_ulim = Inf, p1_xlim = c(1e2, 1e5), p1_ylim = c(1e-2, 1e0),
    d2_ulim = Inf, p2_xlim = c(1e-1, 1e2), p2_ylim = c(1e-2, 1e0),
    d3_ulim = Inf, p3_xlim = c(1e-2, 5e1), p3_ylim = c(1e-2, 1e0),
    p4_xlim = c(1e-3, 5e0), p4_ylim = c(1e-2, 1e0),
)
save_plot("l_j_r_fit_xscale_log", width = 15, height = 10)
"""

l_j_r_fit

l_j_r_der

l_j_r_fit_xscale_log

Plot the median of discontinuity properties with radial distance

Code
function q25(x)
    quantile(x, 0.25)
end

function q75(x)
    quantile(x, 0.75)
end

function stat_info(df; group_cols = [:r, :label])
    cols = [:L_k, :L_k_norm, :j0_k, :j0_k_norm]
    funcs = [median mean std q25 q75]
    
    @chain df begin
        groupby(group_cols)
        combine(cols .=> funcs)
    end
end
stat_info (generic function with 1 method)
Code
### Plot the median of discontinuity properties with radial distance
function plot_median_r(df; add_error=false)
    df_m = stat_info(df)

    fig = Figure()
    plt = data(df_m) * mapping(:r => r_lab, color=:label)

    mappings = [
        ("L_k", l_lab),
        ("L_k_norm", l_norm_lab),
        ("j0_k", j_lab),
        ("j0_k_norm", j_norm_lab)
    ]

    function plt_map(plt, col, label)
        plt *= mapping("$(col)_median" => label)
        lines_map = (visual(Lines) + visual(Scatter))
        errorbars_map = mapping("$(col)_q25", "$(col)_q75") * visual(Errorbars, alpha=0.1, whiskerwidth=10)

        maps = [lines_map]
        add_error && push!(maps, errorbars_map)
        
        plt * reduce(+, maps)
    end

    plts = [plt_map(plt, col, label) for (col, label) in mappings]
    popositions = [(1, 1), (1, 2), (2, 1), (2, 2)] # Define grid positions
    grids = [draw!(fig[pos...], plt) for (pos, plt) in zip(popositions, plts)]
    pretty_legend!(fig, grids[1])

    fig
end
plot_median_r (generic function with 3 methods)
Code
plot_median_r(j_events)
Code
plot_median_r(j_events_taus)
Code
plot_median_r(j_events)
Code
plot_median_r(j_events_taus)
Code
plot_median_r(j_events; add_error=true)

Plot different radial distances in the same plot

Code
function plot_l_r_one(df)
    fig = Figure(size=(1000, 500))

    plt = data(df) * mapping(col=:label, color=:r)
    
    grid1 = draw!( fig[1, 1:2], plt * mapping(:L_k_norm) * density(datalimits=((0, 50),)) )
    grid2 = draw!( fig[2, 1:2], plt * mapping(:L_k) * density(datalimits=((0, 5000),)) )
    legend!(fig[0, 1:end], grid1, titleposition=:left, orientation=:horizontal)

    fig
end

plot_l_r_one(j_events)
Figure 2

Plasma properties with radial distance

Code
plt = data(j_events) * mapping(:radial_distance, :Alfven_speed, color=:label) * (visual(Scatter) + smooth())
plt |> draw
Code
plt = data(j_events) * mapping(:radial_distance, jA_map, color=:label) * (visual(Scatter) + smooth())
plt |> draw
Code
using Statistics
Code
# groupby r and describe the data for each group 
# j_events |> @groupby(_.r) |> @map({r=key(_), j0_k=describe(_.j0_k), L_k=describe(_.L_k)})
@chain j_events begin
    groupby(:r)
    combine(:plasma_density =>  mean, :ion_inertial_length => mean, :b_mag => mean) 
end
5×4 DataFrame
Row r plasma_density_mean ion_inertial_length_mean b_mag_mean
String Float64 Float64 Float64
1 1.0 3.10512 152.217 3.75957
2 2.0 1.42267 242.722 2.10041
3 3.0 0.674452 317.706 1.45649
4 4.0 0.417696 434.665 1.3097
5 5.0 0.343857 569.222 1.27077
Code
plt = data(j_events) * mapping(:radial_distance, :plasma_density, color=:label) * (visual(Scatter) + smooth())
plt |> draw
Code
plt = data(j_events) * mapping(:radial_distance, di_map, color=:label) * (visual(Scatter) + smooth())
plt |> draw

Check discontinuities properties in relation to the local plasma properties

Code
begin
    fig = Figure(size=(1000, 800))

    datalayer = mapping() * visual(Scatter, markersize=1, color=:white, alpha=0.1)

    begin
        layer = histogram(bins=range(1, 5, length=64), normalization=:pdf)
        plt = layer * mapping(di_log_map, l_log_map)
        plt1 = data(j_events) * mapping(layout=:r) * plt
        plt2 = data(w_events) * plt
    
        l_log_limit = ((1.5, 3.5), (1.5, 4.5))
        axis = (;limits=l_log_limit)
        draw!(fig[1:2,1:3], plt1, axis=axis)
        draw!(fig[1:2,4:5], plt2, axis=axis) 
    end

    # Current Density Panels
    begin
        layer = histogram(bins=range(-2, 3, length=64), normalization=:pdf)
        plt = layer * mapping(jA_log_map, j_log_map)
        plt3 = data(j_events) * mapping(layout=:r) * plt
        plt4 = data(w_events) * plt
    
        j_log_limit = ((-1, 3), (-2, 2))
        axis = (;limits=j_log_limit)
        draw!(fig[3:4,1:3], plt3, axis=axis)
        draw!(fig[3:4,4:5], plt4, axis=axis)
    end

    # Make ablines across different r panels
    # begin
        # ablines_data = (; intercepts=[-3,-1,1], slopes=[1,1,1]) 
        # lines = data(ablines_data) * mapping(:intercepts, :slopes) * visual(ABLines, linestyle=:dash)
        # draw!(fig[3:4,1:3], lines)
    # end

    Label(fig[0,1:3], "Juno in Different Radial Distances", fontsize=20)
    Label(fig[0,4:5], "Wind 11 Hz", fontsize=20)

    fig
end
Figure 3
Code
begin
    fig = Figure(size=(1000, 800))

    pdflayer = density() * visual(Contour)
    # small scatter points
    datalayer = mapping() * visual(Scatter, markersize=3)

    # layer = pdflayer + datalayer
    layer = datalayer
    plt = layer * mapping(di_map, l_map)
    plt1 = data(j_events) * mapping(layout=:r) * plt
    plt2 = data(w_events) * plt

    l_limit = ((10^1, 10^5), (10^1, 10^5))
    axis = merge(log_axis, (;limits=l_limit))
    draw!(fig[1:2,1:3], plt1, axis=axis)
    draw!(fig[1:2,4:5], plt2, axis=axis)

    plt = layer * mapping(jA_map, j_map)
    plt3 = data(j_events) * mapping(layout=:r) * plt
    plt4 = data(w_events) * plt

    j_limit = ((10^-1, 10^3), (10^-2, 10^2))
    axis = merge(log_axis, (;limits=j_limit))
    draw!(fig[3:4,1:3], plt3, axis=axis)
    draw!(fig[3:4,4:5], plt4, axis=axis)


    Label(fig[0,1:3], "Juno in Different Radial Distances", fontsize=20)
    Label(fig[0,4:5], "Wind 11 Hz", fontsize=20)

    fig
end
Figure 4
Code
begin
    fig = Figure(size=(1000, 800))
    datalimits_f = x -> quantile(x, [0.05, 0.95])

    begin
        layer = histogram(bins=24, datalimits=datalimits_f, normalization=:pdf)
        plt = layer * mapping(l_norm_log_map, j_norm_log_map)
        plt1 = data(j_events) * mapping(layout=:r) * plt
        plt2 = data(w_events) * plt

        axis = (;)
        draw!(fig[1:2, 1:3], plt1, axis=axis)
        draw!(fig[1:2, 4:5], plt2, axis=axis)
    end

    Label(fig[0, 1:3], "Juno in Different Radial Distances", fontsize=20)
    Label(fig[0, 4:5], "Wind 11 Hz", fontsize=20)

    fig
end

fig

Derivative method validation

Code
fig = plot_l_j_r(j_events_der)
fig
Figure 5